////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//-------------------------------- Ebrahim Foulaadvand, 13 Nov 2012 --------------------------------------//
//                                                                                                        //
//  The routine "NeutronDiffusion" solves the Neutron diffusion equation using forward time centered      //
//  space scheme. Dirichlet boundary condition is used.                                                   //                                     //      //
//                                                                                                        //
//                                                                                                        //
////////////////////////////////////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>

using namespace std;

main()
{


double D=1.,C=1.,L=3.5,tau,h,nL=0.,nR=0.;// D=diffusion constant, C=neutron creation rate, L=system length.

//tau=timestep,h=spatial grid.

int N=60,i,n,Time=20000; //N=number of spatial grids, Time=number of time steps, n=timstep number.

vector<double> ncur(N+1,0),nnew(N+1,0),nAv(Time+1,0);

ofstream file1 ("n n=10.plt");    //output file for temperature profile at timestep n=10.
ofstream file2 ("n n=100.plt");   //output file for temperature profile at timestep n=100.
ofstream file3 ("n n=10000.plt"); //output file for temperature profile at timestep n=1000.
ofstream file4 ("n n=15000.plt"); //output file for temperature profile at timestep n=15000.
ofstream file5 ("nAv.plt");       //output file for the average neutron density.


h=L/N;

cout<<"h= "<<h<<endl; //getch();

//tau=0.5*h*h/(2*D);

tau=0.0005;

//cout<<"tau= "<<tau<<endl; //getch();

nnew[0]=nL; nnew[N]=nR; // implementing Dirichlet boundary condition.

for(i=1;i<N;i++){ncur[i]=0;} ncur[0]=nL; ncur[N]=nR; //implementing the initial condition (here delta function).

ncur[(N/2.0)]=1/h;

// ------------------ Forward time centred space algorithm ----------------------------------


for (n=1;n<=Time;n++){

for(i=1;i<N;i++){ nnew[i]=ncur[i] + (D*tau/(h*h))*( ncur[i+1]+ ncur[i-1]-2*ncur[i] ) + C*tau*ncur[i]; }

if (n==10) for(i=0;i<=N;i++){ file1<<-L/2.0+i*h<<"  "<<nnew[i]<<endl; } //storing in the data file.

if (n==100) for(i=0;i<=N;i++){ file2<<-L/2.0+i*h<<"  "<<nnew[i]<<endl; } //storing in the data file.

if (n==10000) for(i=0;i<=N;i++){ file3<<-L/2.0+i*h<<"  "<<nnew[i]<<endl; } //storing in the data file.

if (n==15000) for(i=0;i<=N;i++){ file4<<-L/2.0+i*h<<"  "<<nnew[i]<<endl; } //storing in the data file.


for(i=1;i<N;i++){ ncur[i]=nnew[i]; }

for(i=0;i<=N;i++){ nAv[n]+=nnew[i]; }

nAv[n]/=(N+1);

file5<<n*tau<<"  "<<nAv[n]<<endl; //storing in the data file.

}

cout<<"end." <<endl;

getch();

}